BACKGROUND OF THE INVENTION 




J^jf frield °f the Invention 

jlfr [0001] The present invention relates to a method of gradual deformation of 
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representations or realizations, generated by sequential simulation, of a model of a 
heterogeneous medium which is not limited to a Gaussian stochastic model, based 
on a gradual deformation algorithm of Gaussian stochastic models. 

Description of the Prior Art 

[0002] In French patent application 98/09,018 a method is described which 
gradually deforms a stochastic (Gaussian type or similar) model of a heterogeneous 
medium such as an underground zone, constrained by a set of parameters relative to 
the structure of the medium. This method comprises drawing a number p (p=2. for 
example) of realizations (or representations) independent of the model or of at least 
part of the selected model of the medium from all the possible realizations and one or 
more iterative stages of gradual deformation of the model by performing one or more 
successive linear combinations of p independent initial realizations and then 
composite realizations are successively obtained possibly with new draws, etc., the 
coefficients of this combination being such that the sum of their squares is 1 . 
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[0003] Gaussian or similar models are well-suited for modelling continuous 
quantity fields but are therefore ill-suited for modelling zones crossed by fracture 
networks or channel systems for example. 

[0004] The most commonly used geostatistical simulation algorithms are those 
referred to as sequential simulation algorithms. Although sequential simulation 
algorithms are particularly well-suited for simulation of Gaussian models, these 
algorithms are not in principle limited to this type of model. 

[0005] A geostatistical representation of an underground zone is formed for 
example by subdividing thereof by a network with N meshes and by determining a 
random vector with N dimensions Z = (Zi, Z2,...Z N ) best corresponding to 
measurements or observations obtained on the zone. As shown for example by 
Johnson, M.E.; in "Multivariate Statistical Simulation" ; Wiley & Sons, New York, 
1987, this approach reduces the problem of the creation of an N-dimensional vector 
to a series of N one-dimensional problems. Such a random vector is neither 
necessarily multi-Gaussian nor stationary. Sequential simulation of Z first involves the 
definition of an order according to which the N elements (Zi, Z2,...Zn) of vector Z are 
generated one after the other. Apart from any particular case, it is assumed that the 
N elements of Z are generated in sequence from Zi to Z N . To determine a value of 
each element Z», (i = 1, .... N), the following operations have to be carried out: 
a) building the distribution of Z|, conditioned by (Zi, Z 2 ..Zm) 
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F c (ZO = P(Z f £*/Zi f Z 2> ...ZM);and 

b) determining a value of Zi from distribution F c (Z ( ). 

[0006] In geostatistical practice, sequential simulation is frequently used to 
generate multi-Gaussian vectors and non-Gaussian indicator vectors. The main 
function of sequential simulation is to determine conditional distributions F c (3) (i = 
1 N), Algorithms and softwares for estimating these distributions are for example 
described in: 

- Deutsch, C.V. et al, "GSLIB (Geostatistical Software Library) and User's Guide" ; 
Oxford University Press, New York, Oxford 1992. 

[0007] Concerning determining the values from distribution F c (3), there also is a 
wide set of known algorithms. 

[0008] The inverse distribution method is considered by means of which a 
realization of Z\\ \ = F c " 1 (Ui) is obtained, where U1 is taken from a uniform distribution 
between 0 and 1 . A realization of vector Z therefore corresponds to a realization of 

vector U whose elements Ui, U 2 U N , are mutually independent and evenly 

distributed between 0 and 1. 
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[0009] It can be seen that a sequential simulation is an operation S which 
converts a uniform vector U = (Ui, U2, ... U N ) to a structured vector Z = (Z 1( Z2,..., 
Z N ): 

Z=S(U) (1). 

[001 0] The problem of the constraint of a vector Z to various types of data can be 

solved by constraining conditional distributions 

F c (Zi) (i = 1 N) and/or uniform vector U = (Ui, U 2 , Un). 

[0011] Recent work on the sequential algorithm was focused on improving the 
estimation of conditional distributions F c (Zj) by geologic data and seismic data 
integration. An article by Zhu, H. et al: "Formatting and Integrating Soft Data : 
Stochastic Imaging via the Markov-Bayes Algorithm" in Soares, A., Ed. Geostatistics 
Troia 92, vol. I : Kiuwer Acad. Publ., Dordrecht, The Netherlands, pp.1 -12, 1993 is an 
example. 

[0012] However, this approach cannot be extended to integration of non-linear 
data such as pressures from well tests and production records, unless a severe 
linearization is imposed. Furthermore, since any combination of uniform vectors U 
does not give a uniform vector, the method for gradual deformation of a stochastic 
model developed in the aforementioned patent application cannot be directly applied 
within the scope of the sequential technique described above. 
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[0013] The method according to the invention thus allows making the two 
approaches compatible, that is to extend the formalism developed in the 
aforementioned patent application to gradual deformation of realizations, generated 
by sequential simulation, of a model which is not limited to a Gaussian stochastic 
model. 

SUMMARY OF THE INVENTION 

[0014] The method allows gradual deformation of a representation or realization, 
generated by sequential simulation, of a model which is not limited to a Gaussian 
stochastic model of a physical quantity z in a heterogeneous medium such as an 
underground zone, in order to constrain the model to a set of data collected in the 
medium by previous measurements and observations relative to the state or the 
structure thereof. 

[0015] The method of the invention comprises applying an algorithm of gradual 
deformation of a stochastic model to a Gaussian vector (Y) having a number N of 
mutually independent variables that is connected to a uniform vector (U) with N 
mutually independent uniform variables by a Gaussian distribution function (G), so as 
to define a chain of realizations u(t) of vector (U), and using these realizations u(t) to 
generate realizations z(t) of this physical quantity that are adjusted in relation to the 
(non- linear) data. 
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[0016] According to a first embodiment, the chain of r alizations u(t) 
of uniform vector (U) is defined from a linear combination of realizations of Gaussian 
vector (Y) whose combination coefficients are such that the sum of their squares is 
one. 

[0017] According to another embodiment, gradual deformation of a number n of 
parts of the model representative of the heterogeneous model is performed while 
preserving the continuity between these n parts of the model by subdividing uniform 
vector (U) into a number n of mutually independent subvectors. 

[0018] The method of the invention finds applications in modeling underground 
zones which generates representations showing how a certain physical quantity is 
distributed in an underground zone (permeability z for example) and which is 
compatible in the best manner with observed or measured data: geologic data, 
seismic records, measurements obtained in wells, notably measurements of the 
variation with time of the pressure and of the flow rate of fluids from a reservoir, etc. 

BRIEF DESCRIPTION OF THE DRAWINGS 

[0019] Other features and advantages of the method according to the invention 
will be clear from reading the description hereafter of a non limitative example, with 
reference to the accompanying drawings wherein: 
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[0020] Fig. 1 shows the medial layer of a realization of a facies model used as a 
reference, generated by sequential simulation of indicatrices, 

[0021] Fig. 2 shows the variation with time of the pressure obtained in a well test 
for the reference model, 

[0022] Figs. 3A to 3E respectively show five initial realizations of the medial layer 
of a reservoir zone, constrained only by the facies along the well, 
[0023] Figs. 4A to 4E respectively show, for these five realizations, the bottomhole 
pressure curves in the reference model compared with those corresponding to the 
initial models, 

[0024] Figs. 5A to 5E respectively show five realizations of the medial layer of the 
facies model conditioned to the facies along the well and adjusted in relation to the 
pressure curve obtained by well tests, 

[0025] Figs. 6A to 6E respectively show, for the five realizations, the bottomhole 
pressure curves in the reference model compared with those corresponding to the 
adjusted models, 

[0026] Figs. 7A to 7E respectively show how the objective functions respectively 
corresponding to these five examples vary with the number of iterations, 
[0027] Figs. 8A to 8E show the gradual deformations generated by an anisotropy 
coefficient change on a three-facies model generated by sequential simulation of 
indicatrices, and 

[0028] Figs. 9A to 9E show the local gradual deformations of a three-facies model, 
generated by sequential simulation of indicatrices. 
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DETAILED DESCRIPTION OF THE INVENTION 



[0029] A study zone is considered that is subdivided by an N-mesh grid. 
Realizations or representations of a stochastic model of a certain physical quantity z 
representing for example the permeability of the formations in the zone are attempted 
to be made. The model that is sought must adjust to data obtained by measurements 
or observations at a certain 

[0030] Adjustment of a stochastic model to non-linear data can be expressed as 
an optimization problem. The quantity f 0 * = (fi obs f f 2 0b3 f f 3 obs ... f P obs ) designates the 
vector of the nonlinear data observed or measured in the studied medium (the 
reservoir zone), and the quantity f = (fi, f 2 , h .. ■ fp) is the corresponding vector of the 
responses of the stochastic model of the permeability Z = (Zi, Z 2) Z N ,). The 
problem of constraining the stochastic model of Z by observations consists in 
generating a realization z of Z which reduces to a rather low value an objective 
function that is defined as the sum of the weighted rms errors of the responses of the 
model in relation to the observations or measurements in the reservoir zone, i.e.: 
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wher cuu r presents th w ight assigned to respons fj. Functions f| (i=1, 2, p) 
and objective function O are functions of vector Z. This presents an optimization 
problem of dimension N. 

[0031] In order to extend the formalism developed in the aforementioned patent 
application to the gradual deformation of realizations generated by, but not 
necessarily limited to Gaussian sequential simulation, a starting point is from a 
Gaussian vector with N variables Yi, with i = 1, 2, N, mutually independent, of 
zero mean and of variance equal to 1, and N mutually independent uniform variables 
Ui, U 2 , U 3 ,...U N are defined by: 

U, = G(Y,) Vi = 1,2, „.,N 

where G represents the standardized Gaussian distribution function. 

[0032] Assuming this to be the case, the gradual deformation algorithm developed 
within a Gaussian frame is applied to the Gaussian vector Y = (Yi, Y2,..., Y N ) in order 
to build a continuous chain of realizations of uniform vector 
U (Ui, U 2 ,... f U N ). Given two independent realizations y A and y b of Y, the chain of 
realizations u(t) of vector U obtained with the following relation is defined: 

u(t) = G(y a , cos t + y b , sin t) (2). 
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[0033] For each t, u(t) is a realization of vector U. A vector z(t) which is, for each t, 
a realization of random vector Z is then obtained by sampling of the conditional 

distribution F c (Zi) (i=1, 2 N) using the elements of vector u(t). Parameter t can 

consequently be adjusted as in the Gaussian case so as to adjust z(t) to non-linear 
data. This procedure is iterated until satisfactory adjustment is obtained. 

Adjustment of a Facies Model to Pressure Data Obtained by Means of Well 
Tests 

[0034] In order to illustrate application of the stochastic optimization method 
defined above, adjusting a stochastic reservoir model to pressure data obtained by 
means of well tests is attempted. Building of the reservoir model is derived from a 
real oil formation comprising three types of facies: two reservoir facies of good quality 
(facies I and 2) and a reservoir facies of very bad quality (facies 3). Table I defines 
the petrophysical properties of the three facies: 





K x (md) 


K v (md) 


K z (md) 


<P (%) 


c, (1 0 s bar 1 ) 


Facies 1 


10 


10 


10 


17 


2. 1857 


Facies 2 


1 


1 


1 


14 


2.0003 


Facies 3 


0.1 


0.1 


0.001 


9 


1.8148 



[0035] In order to represent the specific facies distribution of the oil formation, a 
binary realization is first generated to represent facies 3 and its complement. Then, in 
the complementary part of facies 3, another binary realization independent of the first 
one is generated to represent facies 1 and 2. The formation is discretized by means 
of a regular grid pattern of 60x59x15 blocks 15rnx15mx1.5m in size. An exponential 
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variogram model is used to estimate the conditional distributions. The main 
anisotropy direction is diagonal in relation to the grid pattern. The ranges of the 
variogram of facies 3 in the three anisotropy directions are 300m, 80m and 3m 
respectively. The ranges of the variogram of facies 1 and 2 in the three anisotropy 
directions are 150m f 40m and 1. 5m respectively. The proportions of facies 1, 2, 3 are 
6%, 16% and 78% respectively. 

[0036] A well test has been carried out by means of a finite-difference well test 
simulator as described by: Blanc, G. et al: "Building Geostatistical Models 
Constrained by Dynamic Data - A Posteriori Constraints" in SPE 35478, Proc. 
NPF/SPE European 3D Reservoir Modelling Conference, Stavenger, Norway, 1996, 

[0037] The medial layer of a realization used as the reference model for this 
validation can be seen in Fig. 1. The section of the well that has been drilled runs 
horizontally through the medial layer of the reservoir model along axis x. The 
diameter of the well is 7.85cm, the capacity of the well is zero and the skin factors of 
facies 1 , 2 and 3 are 0, 3 and 50 respectively. The synthetic well test lasts for 240 
days with a constant flow rate of 5 m 3 /day so as to investigate nearly the entire oil 
field. Fig. 2 shows the pressure variation with time. 

[0038] The objective was to build realizations of the facies model constrained by 
the facies encountered along the well and by the pressure curve obtained during well 
testing. The objective function is defined as the sum of the rms differences between 
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the pressure responses of the reference model and the pr ssure responses of the 
realization. Since the dynamic behavior of the reservoir model is mainly controlled by 
the contrast between the reservoir facies of good and bad quality, the binary 
realization used to generate facies 1 and 2 has been fixed first and only the binary 
realization used to generate facies 3 has been deformed for pressure data 
adjustment. 

[0039] The pressure responses resulting from the well tests for the five 
realizations of Figs.3A to 3E are different from that of the reference model, as shown 
in Figs. 4A to 4E. Starting respectively from these 5 independent realizations, by 
using the iterative adjustment method above, after several iterations, five adjusted 
are obtained realizations (Figs. 5A to SE) for which the corresponding pressure 
curves are totally in accordance with those of the reference model, as shown in Figs. 
6A to6E. 

Gradual Deformation in Relation to the Structural Parameters 

[0040] In many cases, sufficient data for deducing the structural parameters of the 
stochastic model: mean, variance, covariance function, etc., are not available. These 
structural parameters are often given in terms of a priori intervals or distributions. If 
their values are wrong, it is useless to seek a realization adjusted to non-linear data, 
It is therefore essential for applications to be able to perform a gradual deformation of 
a realization with simultaneous modification of random numbers and structural 
parameters. The sequential simulation algorithm defined by equation (1) makes 
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possible changing, simultaneously or s parat ly, a structural operator S and a 
uniform vector U. Figs. 8A to 8E show the gradual deformations obtained for a fixed 
realization of uniform vector U when the anisotropy coefficient is changed. 

Local or Regionalized Gradual Deformation 
[0041] When the observations are spread out over different zones of a studied 
formation, an adjustment using global deformation would be ineffective because the 
improvement obtained in a zone could deteriorate the improvement in another zone. 
It is therefore preferable to apply gradual deformations zone by zone. Consider a 
subdivision of vector U into a certain number n of mutually independent subvectors 
U 1 , U 2 ,..., U n , which allows performing their gradual deformation individually. 
Separate application of the gradual deformation algorithm to each subvector U\ U 2 ..., 
U n allows obtaining a function of dimension n of uniform vector U: 









'G(Xl cos r^'sin/,) " 




u\t 2 ) 




G(Y a 2 cos t 2 + Y? sin t 2 ) 








_G(Y a " cos t„ + Y? sin t„)_ 



where Y'a and Yb for any i = 1 , 2,..., n, are independent Gaussian subvectors. For a 
set of realizations of Y* a , and Yb, a problem of optimization of n parameters ti, t 2f ..., tn 
is solved to obtain a realization that improves or maintains the data adjustment. This 
procedure can be iterated until satisfactory adjustment is obtained. 
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[0042] Gradual local deformations thus allow significant improvement of the 
adjustment speed in all the cases where measurements or observations are spread 
out over different zones of the medium. 

[0043] The effect of this gradual local deformation on the three-facies model of 
Figs. 9A to 9E can be clearly seen therein where only the enclosed left lower part is 
affected. 

[0044] The method according to the invention can be readily generalized to 
gradual deformation of a representation or realization of any stochastic model since 
generation of a realization of such a stochastic model always comes down to 
generation of uniform numbers. 
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